NASA/CR- 1998-208963 
ICASE Report No. 98-55 



A Gas-kinetic Scheme for Reactive Flows 


Yongsheng Lian 

The Hong Kong University of Science and Technology, Hong Kong 
Kun Xu 

ICASE, Hampton, Virginia 
and 

The Hong Kong University of Science and Technology, Hong Kong 


Institute for Computer Applications in Science and Engineering 
NASA Langley Research Center 
Hampton, VA 


Operated by Universities Space Research Association 



National Aeronautics and 
Space Administration 

Prepared for Langley Research Center 
under Contract NAS 1 -97046 


Langley Research Center 
Hampton, Virginia 23681-2199 


December 1998 



Available from the following: 


NASA Center for AeroSpace Information (CASI) 
7121 Standard Drive 
Hanover, MD 21076-1320 
(301)621-0390 


National Technical Information Service (NTIS) 
5285 Port Royal Road 
Springfield, V A 22161-2171 
(703) 487-4650 


1 1 r 



A GAS-KINETIC SCHEME FOR REACTIVE FLOWS * 

YONGSHENG LTAN t AND KUN XU * 


Abstract. In this paper, the gas-kinetic BGK scheme for the compressible flow equations is extended to 
chemical reactive flow. The mass fraction of the unburnt gas is implemented into the gas kinetic equation 
by assigning a new internal degree of freedom to the particle distribution function. The new variable can be 
also used to describe fluid trajectory for the nonreactive flows. Due to the kinetic governing equation, the 
current scheme basically solves the Navi er- Stokes chemical reactive flow equations. Numerical tests validate 
the accuracy and robustness of the current kinetic method. 
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Subject classification. Applied Numerical Methods 

1. Introduction. There are mainly two numerical approaches to the solution of the compressible Euler 
equations, namely, the Godunov and the Boltzmann schemes. Broadly speaking, the Godunov scheme is 
based on the Riemann solution and characteristics play an important role in the description of the gas 
evolution. However, the Boltzmann scheme uses the microscopic particle distribution function as the basis 
in the construction of the flux function and the Euler solution is considered as a limiting case when the 
particle collision time goes to zero. The Godunov and the Boltzmann schemes are based on two different 
physical interpretations of flow motion. Due to the possible implementation of nonequilibrium gas property 
in the kinetic scheme, both the robustness and accuracy of the scheme can be maintained [10]. 

In this paper, we extend the gas-kinetic BGK scheme for the nonrcactive compressible Euler equations 
to the reactive flows. In order to implement the mass fraction into the kinetic formulation, one new internal 
degree of freedom z is implemented in the gas distribution function. For nonreactive flows, this function can 
be also used for the tracking of fluid trajectory. In the reactive flow calculations, wc are only accounting for 
two species, which are the unburnt and burnt gases. The unburnt gas is converted to burnt gas via a simple 
irreversible reaction process. As a special application, the new scheme is used in the study of detonation 
waves in both 1-D and 2-D cases. 

The inviscid reacting compressible Euler equations in 1-D case arc 

Pt + (pU) x — 

(p U) t + (pU 2 + p) x = 0, 

0 oZ)t + (pZU) x = —pK(T)Z, 

(pt)t + (pfU + pU)x = q 0 pK(T)Z , 

where p is the density, U the velocity, p the pressure, Z the mass fraction of unburnt gas, and qo is the 
amount of heat released per unit mass by reaction. The total energy density is pe = \pU 2 4- pe, where pe is 
the internal energy. We assume that both unburnt and burnt gases have the same 7. The equation of state 
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can be expressed as p — pTR/m , where R is the gas constant and m is the molecular mass. K(T) is the 
chemical reactive rate, which is a function of temperature. The specific form of K{T) will be given in the 
numerical section. 

Many researchers have been working on the numerical solution of equations (1.1). A partial list of 
references includes Colella et.al. [2], Lindstrom [6], Engquist and Sjogreen [3], and Jeltsch and Klingcnstein 
[4]. Mostly, a splitting scheme is used to solve the above equations and the flow variables inside each cell are 
updated through 

dW 1 

~dT = ^ {F ^ 2{t) “ F W*)) + W), 

where Wj = (p, pU, pZ, pc)J is the cell-averaged conservative variables, S = (0,0, — pK(T)Z, qopK(T)Z)J is 
the source term, and the flux function is obtained by solving Eq.(l.l) without considering the source 

term. In this paper, a gas-kinetic model and the corresponding kinetic scheme for the evaluation of the flux 
function 1/ 2 of the homogeneous Euler equations will be presented, and the source term S(Wj ) in the 
above equation is treated implicitly for the update of flow variables Wj inside each cell. 

2. A Gas-Kinetic Model. A gas-kinetic BGK model for Eq.(l.l) without the source terms can be 
constructed as the following equation, 

(2.1) ft + uf x = 

where / is the gas-distribution function, u the particle velocity, and Q(f, f) = (g — f)/r the particle collision 
term [1]. The equilibrium state g has the form, 

K + 2 

g = p(^j e-Hiv-ufHt-zf+e), 

where K is the number of dimensions of the internal variable £ and is related to 7, 


K = { 3-7)/(7 -l), 


and £ 2 = £ 2 + £ 2 + ... + Ck- The term A is a function of the gas temperature T with the relation A = m/ 2 KT 
and k is the Boltzmann constant. 

The connection between the distribution function / and the macroscopic flow variables is 

(p, pU, pZ, pe) T = / ip a fdudzd£, 

where d£ = d^\d^2--d^K and 

tp a = (1 ,u,z, i(u 2 + 4 2 )) T 

are the moments for density p, momentum pU> mass fraction pZ, and total energy pc. The fluxes for the 
corresponding macroscopic variables are 

(Fpy FpU , FpZ , F pe ) T = / uifcafdudzd 

For the homogeneous flow equations (1.1) without the source terms, the compatibility condition of the 
collision term in the Boltzmann equation is 


( 2 . 2 ) 


/ 


Q(/; f)ip a dudzd£ = 


0 

0 

\o/ 



For the equilibrium flow with / = g> the homogeneous Euler equations with the inclusion of mass fraction 
can be recovered by taking the moments of xp a to Eq.(2.1), 


/ 


1 

u 

z 
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(g t -f- ug x )dudzd £ = 0, 


\\{u 2 +e)J 


and the resulting equations become 
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So, the corresponding pressure is p = p/2 A and the internal energy density pe goes to 

(K + l)p 


pe = 


4 A 


To the first order of r, the Chapman-Enskog expansion gives 


/ = 9 - r{g t + ug x ), 


and the BGK model automatically reduces to the Navier-Stokes equations, 


p A 


( pu \ 
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where the dynamic viscous coefficient is 77 = rp. In the 2-D cases, similar viscous governing equations can be 
derived from the BGK model [ 10 ]. As a result, for the chemical reactive flows, the real governing equations 
solved by the kinetic BGK scheme are the reactive Navier-Stokes equations instead of inviscid equations 
(1.1). This is basically one of the direct reason for the robustness of kinetic BGK scheme. 

Remark: in the above equations, the function Z has no dynamical effect on the gas evolution, it only provides 
additional information about the flow property, which can be the mass fraction for the reactive flow, level 
set function for the interface tracking, color function for the fluid trajectory capturing, and the pollution 
concentration for certain gas species. 


3 . Gas-Kinetic Flow Solver . In order to evaluate the numerical fluxes across a cell interface #7+ 1/2, 
we need to get the gas distribution function there. The general solution of / at the cell interface #7+1/2 and 
time t is 

1 ft / 

( 3 . 1 ) f(x j+ i/ 2 ,t,u,z ,0 = - / g{x',t’ ) u,z^)e^ (t ~ t ' )/T dt' + e~ t/r f 0 {x j+1/2 - ut), 

r Jo 

where x f = #7+1/2 — u(t — t f ) is the trajectory of the particle motion and /o is the initial gas distribution 
function / at the beginning of each time step (t — 0 ). Two unknowns, g and /o in Eq.( 3 . 1 ), have to be 
addressed in the above equation in order to obtain the explicit form of /. 

Generally, the distributions of /o and g around the cell interface #7+1/2 and time t ~ 0 are obtained 
using the Taylor expansion of the Maxwellian distribution function, for example 


( 3 . 2 ) 


f 0 


g l (l +a l {x - x j+1/2 )) , x<Xj +\/2 

g r (l+a r (x -x j+1/2 )), x > Xj +l/2 
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and 


(3.3) g = so (1 + (1 — H[x - x j+1 /2 \)a l {x - x j+1/2 ) + H[z - x j+l/2 ]a r {x - x j+l/2 ) + At ) , 

where g\g r and g Q are local Maxwellians located at the left, right and middle of a cell interface. The 
parameters a l ,a r ,a l ,a r have the following form 

a = oi + a 2 u + a$z + a±(u 2 + £ 2 ), 


and all parameters (ai, 02 , « 3 , < 14 ) can be found from the slopes of the corresponding macroscopic variables. 
H[x] is the Heaviside function defined as 



x < 0 
x > 0 


The reason and detailed formulation in the determination of /o is presented in [12, 10]. The only difference 
here is that we need to use the macroscopic distribution pZ in the determination of the a$z term. 

After /o is obtained, the equilibrium state £0 at a cell interface 

K + 2 

go = Po(-) 2 

7 r 

is determined as follows. Taking the limit of t — > 0 in Eq.(3.1) and substituting its solution into Eq.(2.2), 
the compatibility constraint at (x = Xj+ if 2 , t = 0) gives 


Wn = 


/ goip a dudzd€ = / / g l ^ a dudzd^ + / / g r ip a dudzd £. 

J u>0 J J ti<0 J 


Similarly the corresponding slopes of g in Eq.(3.3) can be obtained from the macroscopic slopes between the 
cell averaged flow quantities Wj and W J+ 1 and the above value Wq at the cell interface [12]. 

After substituting Eq.(3.2) and Eq.(3.3) into Eq.(3.1), the final gas distribution function at a cell interface 
is 


f{xj+i/ 2 ,t,u,z,£) = (1 - e t/T )g 0 + (t{- 1 + e t/r )+te t/T ) (a*H[u] + a r (l - H[u])) ug 0 
(3.4) +r(t/T - 1 + e~ t/T )Ag 0 

+e~ t / T ((1 — uta l )H[u]g l + (1 — uta r )( 1 — H[w])p r ) . 


The only unknown in the above equation is A term, which is determined by implementing the compatibility 
condition over the whole time step At at the location x J+1 / 2 , 


rAt 

Jo 


t, u, 2 , £) — /(x J+ i/ 2 i u, z , £))dtdudzd£ = 0. 


There is no iteration involved in the determination of A from the above equation [12]. After / is obtained, 
the time- dependent numerical fluxes in the x-direction across the cell interface can be computed as 


(3.5) 


/*> \ 
f P u 
F P z 

V Fpt ) j+1/2 


/- 
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\ 


f(z j+l/2 ,t,u,z, £)dudzd£. 


\u^+e)/ 


By integrating the above equation over the whole time step At, we get the total W = (p, pU , pZ, pe) transport. 


IT! 
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4. Numerical Examples. In the numerical examples reported in this section, the van Leer limiter is 
used for the reconstruction of conservative variables W at the beginning of each time step. Unless specifically 
stated, the gas constant 7 is equal to 1.4. The first three test cases are about the nonreactive flow and the 
mass fraction is used as an interface tracer for the fluid evolution; the following two test cases are about 1-D 
and 2-D detonation wave calculations. 

Case(l): Diffusion of mass fraction function 

As analyzed in section 2, the real governing equation obtained from the kinetic BGK model for the 
function pZ is the advection diffusion equation, 

(pZ)t + (pUZ) x = T(^Z x ) x . 

In order to test the above governing equation, we set two uniform initial flow conditions with 
(4.1) (p = l,p = 1,C/ = 0), and (p = l,p = 1, U = 0.5). 

The computational domain consists of 200 grid points with cell size Ax = 1.0. The function Z is initially 
assigned with the value 

z r-i. *<ioo, 

\l. x > 100. 

Two fixed collision times r = 0.03 and r = 0.015 are used in the computations, which correspond to 
viscosity coefficients u = 0.03 and v = 0.015 respectively. At the output time t = 100, the numerical and 
exact solutions for both cases are shown in Fig, (5.1). The results confirm that the BGK scheme does solve 
the advection diffusion equation for the mass fraction function Z. 

Case(2) Fluid trajectory in the shock tube case [9] 

The forward-facing case is carried out on a uniform mesh of 120 x 40 cells and Ax = Ay = 1/40. We 
choose the color function Z at the inlet x = 0 with the following boundary conditions 

Z= < 


The computed density and pressure distributions are presented in Fig. (5.2). In the same figure, the contours 
of function Z arc added, from which the interfaces between different “colored” fluid and the fluid trajectories 
can be clearly observed. For example, the fluid particles change direction after passing through the oblique 
shock. 

Case(3) Rayleigh- Taylor instability [7, 5] 

This computation is performed on a rectangular domain of x € [0, 1] and y e [0, 2] with reflecting 
boundary conditions on the lower and upper sides of the domain and periodic ones in the horizontal direction. 
The gravity is directed downward with dimensionless gravitational constant G = 0.5. 

The densities next to the initial fluid interface at y = 1 are p\ — 0.5 and p 2 = 1.0 with the ratio 
P2/P1 =2:1, and the functions Z are 1.0 below the interface and —1.0 above that. The value of the pressure 
at the fluid interface (y = 1) is 1/1.4, and isothermal conditions are used to determine flow distributions in 
both the upper and lower parts. The initial density perturbation at the interface is added with the form 
Sp — 0.05(1 — cos(27rx)). Since the heavy fluid is located on top of the light fluid, it stays in an unstable 


1.0 for0 <y<^, 

- 1.0 for ^ < y < 

1-0 for § < 2/ < I, 
- 1.0 for H < y < 1 . 
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situation when the system is subjected to gravity. The computed contours of function Z with the values 
Z = [“0.5, 0.0, 0.5] at output time t = 10.0 on three different mesh sizes (Ax = Ay — 1/64, 1/128, 1/192) arc 
shown in Fig.(5.3). Since the collision time r = 4 x 10 -4 is fixed in all these three cases, the physical viscosity 
coefficient rpj p keeps the same value. Therefore, even with the mesh refinement the simulation results are 
basically identical. If the Riemann solution of the inviscid Euler equations is used in the flux evaluation 
1 / 2 , the numerical results usually do not converge with the mesh- refinement [7]. The nonconvergence 
of the numerical results is more serious for the detonative reactive flows using the exact Godunov method 
[ 8 , 6 ]. 

Case(4) 1-D denotation wave 

In this case, we are going to study the formation of the ZND wave for the following reaction kinetics, 


(4.2) 


K(T) = 


■Ko> T > To, 
0, T < T 0 , 


where T 0 is the ignition temperature and K 0 the reaction rate. This specific case is taken from [4]. The 
initial data are piece- wise constant, which defines the Chapman- Jouget detonation wave: 

(pi = 1.9690 x lO" 3 -^, Ui = 4.8057 x 10 4 — , pi = 7.9434 x 10 6 — Z = 0), 
cm 3 sec cmsec 1 

(, Pr = 1-2010 x lCT 3 ^, Ui = 0.0, pi = 0.8321 x 10 6 ^^, Z = 1.0). 


In the calculation, the gas constant R is equal to 8.3143 x 10 7 cm 2 g j sec 2 K mol, the molecular weight m = 
36 gjmol, the ignition temperature To = 500°1T, the reaction rate Kq = 0.582458 x 10 lo /s, and the heat 
release q 0 = 6.9283 x I0 9 cm 2 /sec 2 . The spatial step size Ax used in each case is varied according to 
Ax = cxRqi where Rq = 5.347 x 10~ 6 cm and the parameter a takes the values 0.01, 0.1, 1 in the three cases. 
The results at subsequent times with different a arc shown in Fig.(5.4)-(5.6). From these figures, wc find 
that a detonation wave is a strong shock wave propagating into a reactant, followed by a thin zone of reaction 
which supports the shock. 

Case(5) 2-D denotation wave [6] 

As illustrated in [6], with the mesh refinement the reactive Euler solvers can generate unphysical solutions 
in the complicated oscillating detonation waves. The spurious solution appears even using the exact Godunov 
method [8]. 

The initial condition for the 2D denotation simulation is an exact traveling solution of the ZND wave. 
The reaction rate K(T ) has the following Arrhenius formulation, 


(4.3) K(T) = K Q T a e~ E ' T . 

The parameters used are qo = 50, E = 50,7 = 1-2 and a = 0. The reaction rate Kq is set to be 10 4 . 
The initial data is a one- dimensional ZND profile in the x-direction. The ZND wave connects the left state 
pi = 1.731379, Ui = 3.015113 Vj = 0, piei = 130.4736, Zi = 0 by a Chapman-Jouget detonation with the 
right state p r = 1, U r = 0, V R = 0, p r e r = 15, Z r = 1. The computational domain is 0.6 x 1.0, and the cell 
size used is Ax = Ay = 1/400. A periodic perturbation is imposed in the y-direction of the initial ZND 
profile, where the initial data W{x , y, 0) is set to W Z nd{x + AxNINT(^cos(47ry))), where NINT(z) is the 
nearest integer close to z. The simulation results around the ZND wave front for the subsequent times from 
5/64,6/64, ..., 16/65 are shown in the Fig.(5.7). From this figure, we can clearly see the oscillating profile of 
the ZND wave front and the “explosion within explosion ” phenomena due to the collision of triple points. 

For the Godunov scheme, due to the inadequate dissipation in the gas evolution stage, the shock insta- 
bility and carbuncle phenomena are intrinsically rooted [11], and the robustness of the Godunov scheme can 
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hardly be achieved for the complicated flow system in the high resolution calculations, such as the chemical 
reactive and MHD equations. A detail analysis of the dissipative mechanism in the Godunov method is 
presented in [11]. However, for the gas-kinetic BGK scheme, we are basically solving the viscous governing 
equations even for the inviscid target equations. So, the robustness is well maintained. This kind of approach 
is physically founded, because we are obtaining the numerical solutions on the discretized space and time, 
where the spatial and temporal resolution is limited by the cell size and time step. The subcell smearing is 
equivalent to the intrinsic dissipation. In order to remove the unphysical solutions in Riemann solver based 
methods, the reactive Navier-Stokes equations arc solved directly in [6]. However, the evaluation of the vis- 
cous terms requires substantial computation resources. In this aspect, the gas-kinetic scheme of the current 
paper is efficient since the viscous and heat conduction terms have been included in the gas distribution 
function (3.5) already. 

5. Conclusion. In this paper, we have extended the BGK scheme to the chemical reactive flow with 
the inclusion of one more internal degree of freedom in the gas distribution function to account for the mass 
fraction. This mass function can be also used to track the fluid interfaces for the nonreactive flows. The 
numerical results confirm the robustness, accuracy and efficiency of the BGK method. 
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Fig. 5.1. Propagation of function Z with velocity U = 0 (top) and U = 0.5 (bottom). The collision times used are r — 0.03 
(left) and r — 0.015 (right) respectively. The solid lines are exact solutions and the circles are numerical solutions. 




Fig. 5.2. Density and pressure contours for the step case. The contours of function Z with Z — [—1/32, 0, 1/32] are added 
in these plots to show the fluid trajectory. The mesh size used here is 120 x 40. 
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64*128 points 128*256 points 192*384 points 

Fig. 5.3. In the Rayleigh- Taylor instability case, the interface between the heavy and light fluid is captured with the 
help of function Z, and the contours have the values of Z = [-0.5, 0.0, 0.5], From left to right, the mesh sizes used are 
64 x 128, 128 x 256, 192 x 384 and the collision time in each case keeps the same value r = 4 x 10~ 4 . Therefore, the physical 
viscous coefficient is the same , so are the simulation results. 
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Fig. 5.5. Density, pressure , velocity and mass fraction Z plots at time steps 0,2000, 10000 with a — 0.1. 





Fig. 5.6. Density, pressure , velocity and mass fraction Z plots at time steps 0, 2000, 10000 with a — 1. Due to the 
large cell size in this case, the peak values are reduced and the profiles get smeared in comparison with Fig. (5.4) and Fig. (5.5). 
If ot is continuously increasing, spurious solutions of one cell per time step will appear [2]. 




Fig. 5.7. Density distribution of the propagating detonation front at time t = 5/64,6/64, ...,16/64 (from left — ► right, top 
— *■ bottom). The phenomena of “explosion within explosion ” can be clearly observed at the leading shock front 
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